Advanced_Analytics_for_Operational_Planning_in_Bakeries

Vorgehensmodell

Business Understanding

Der Bäckereibetrieb befindet sich im Schwarzwald im Kreis Rottweil und hat insgesamt 13 Filialen in 8 verschiedenen Städten. Während viele Discounter und Bäckereiketten die Waren in ihrer Filiale im Ofen aufbacken, wird bei der Brantner Bäck GmbH & Co. KG alles frisch in der Produktionsstätte auf dem Hardt produziert. Von dort werden dann die einzelnen Filialen beliefert. Für die Lieferung muss vorab entschieden wer- den wie viele Waren die einzelnen Filialen benötigen.

Für die zukünftige Bestimmung der optimalen Produktionsmenge soll ein datengetriebener Ansatz verfolgt werden. Aus diesem Grund sollen nun vermehrt die Daten in den Vor- dergrund gerückt werden, um bessere Entscheidungen zu treffen. Dafür benötigt es die richtige Datengrundlage, sowie aussagekräftige Insights. Deshalb ist es Ziel ein Sys- tem zu entwickeln, das die Bäckerei bei ihren Entscheidungen unterstützt. Insgesamt wurden für das Projekt die folgenden drei Ziele definiert. • Transparenz zur besseren Entscheidungsfindung • Optimale Produktions- und Liefermengen • Komplette Automatisierung Das System sollte die Bäckerei bei der Entscheidungsfindung unterstützen und Trans- parenz in die Daten bringen. D

Data Understanding

Packages laden

library(tidyverse)
library(plotly)
library(readxl)
library(dplyr)
library(ggplot2)
library(tidyverse)
library(readxl)
library(prophet)
library(data.table)
library(dplyr)
library(ggplot2)
library(onehot)
library(mltools)
library(flexdashboard)
library(randomForest)

Daten einlesen

gc()

##           used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 2018146 107.8    3775778 201.7  3137811 167.6
## Vcells 3679622  28.1    8388608  64.0  5586986  42.7

data <- read.csv("~/R 3.5 Files/Masterarbeit/Brantner_Daten_1507.csv")

Key Statistics

#Zeilenanzahl
nrow(data)

## [1] 1264645

#Spaltenanzahl
ncol(data)

## [1] 26

#Übersicht der Features
head(data, 0)

##  [1] X              filiale        datum          artnr         
##  [5] artbez         anzahl         anfangsbestand endbestand    
##  [9] summe          std5           std6           std7          
## [13] std8           std9           std10          std11         
## [17] std12          std13          std14          std15         
## [21] std16          std17          std18          std19         
## [25] std20          id            
## <0 rows> (or 0-length row.names)

#Datum als Datum speichern
data$datum <- as.Date(data$datum)

#Zeitspanne
min(data$datum)

## [1] "2016-10-01"

max(data$datum)

## [1] "2019-07-14"

max(data$datum) - min(data$datum)

## Time difference of 1016 days

#Anzahl an Produkten
length(unique(data$artnr))

## [1] 961

#Anzahl an Filialen
length(unique(data$filiale))

## [1] 16

#3 davon gibt es nicht mehr

#data$artnr <- as.character(data$artnr)
data$filiale <- if_else(data$filiale == "Weilerstr", as.character("Weilerstr."), as.character(data$filiale))

data[is.na(data)] <- 0

Umsatzanalyse - Umsatz nach Wochentagen

#Tag hinzufügen
ag <- aggregate(data$summe, by=list(data$datum), FUN=sum)
ag$x <- as.numeric(ag$x)
ag$Group.1 <- as.Date(ag$Group.1)
colnames(ag) <- c("datum","summe")

ag$day <- weekdays(ag$datum)

ax <- list(
  title = "",
  zeroline = FALSE,
  showline = FALSE,
  showticklabels = FALSE,
  showgrid = FALSE
)

plot_ly(data = ag, x = ~datum, y = ~summe, color = ~day, hoverinfo = "none") %>%
  layout(xaxis = list(showticklabels = FALSE), yaxis = list(showticklabels = FALSE))

Umsatzanalyse - Umsatzverteilung nach Produkten

#Relevante Produkte und Filialen
prods <- c(1525, 1108, 1101, 4129, 1526, 1528, 4136, 38, 2120)
fil <- c("Weilerstr.", "Rath.Platz")

#Filtern
df_plot <- data %>%
  filter(artnr %in% prods) %>%
  filter(filiale %in% fil)

prods_bez <- unique(df_plot$artbez)

#Als Faktor speichern
df_plot$filiale <- as.factor(df_plot$filiale)
df_plot$artbez <- as.factor(df_plot$artbez)
df_plot$artnr <- as.character(df_plot$artnr)
df_plot$artnr <- as.factor(df_plot$artnr)

#Farbpalette erstellen
library(RColorBrewer)
mycolors <- brewer.pal(n = 9, name = "Blues") 
mycolors <- c(mycolors, "#001642", "#000019")
mycolors <- rev(mycolors)
mycolors <- mycolors[1:9]


plot_ly(df_plot, y = ~summe, x = ~filiale, color = ~artnr, type = "box", colors = mycolors, hoverinfo = "none") %>%
  layout(boxmode = "group", legend = list(orientation = 'h'), xaxis = list(showticklabels = FALSE), yaxis = list(showticklabels = FALSE))

Retourenanalyse - Datenvorverarbeitung

Brezeln & Laugenstangen haben jeweils eine Buttervariante auf einer anderen Artikelnummer. Diese müssen zusammengefügt werden.

#Filtern auf releavante Produkte & Filialen
prods <- c(1525, 1108, 1101, 4129, 1526, 1528, 4136, 38, 2120)
fil <- c("Weilerstr.", "Rath.Platz")

#filtern
df_plot <- data %>%
  filter(artnr %in% prods) %>%
  filter(filiale %in% fil)

#Brezeln & Butterbrezeln zusammenfügen
brezeln <- c(1525, 4129)
df_brezeln <- df_plot %>%
  filter(artnr %in% brezeln )

#Artikel gruppieren und aggregieren
df_brezeln <- df_brezeln %>% 
  group_by(datum, filiale) %>% 
  summarise(anzahl = sum(anzahl), anfangsbestand = sum(anfangsbestand), endbestand = sum(endbestand), summe = sum(summe), std5 = sum(std5) , std6 = sum(std6), std7 = sum(std7) , std8 = sum(std8), std9 = sum(std9), std10 = sum(std10), std11 = sum(std11), std12 = sum(std12), std13 = sum(std13), std14 = sum(std14), std15 = sum(std15), std16 = sum(std16), std17 = sum(std17), std18 = sum(std18), std19 = sum(std19), std20 = sum(std20))

#Benennen
df_brezeln$artnr <- 1525
df_brezeln$artbez <- "Brezeln und Butterbrezeln"



#Laugenstangen & Butterlaugenstangen zusammenfügen
laugenstange <- c(1528, 4136)
df_laugenstange <- df_plot %>%
  filter(artnr %in% laugenstange )

#Artikel gruppieren und aggregieren
df_laugenstange <- df_laugenstange %>% 
  group_by(datum, filiale) %>% 
  summarise(anzahl = sum(anzahl), anfangsbestand = sum(anfangsbestand), endbestand = sum(endbestand), summe = sum(summe), std5 = sum(std5) , std6 = sum(std6), std7 = sum(std7) , std8 = sum(std8), std9 = sum(std9), std10 = sum(std10), std11 = sum(std11), std12 = sum(std12), std13 = sum(std13), std14 = sum(std14), std15 = sum(std15), std16 = sum(std16), std17 = sum(std17), std18 = sum(std18), std19 = sum(std19), std20 = sum(std20))

#Benennen
df_laugenstange$artnr <- 1528
df_laugenstange$artbez <- "Laugenstange und Butterlaugenstange"

#Aus dem vorherigen Dataframe selektieren
df_plot <- df_plot %>%
  select(-X,-id)

#Entsprechende Spaltennamen
header <- colnames(df_plot)
df_brezeln <- df_brezeln[header]
df_laugenstange <- df_laugenstange[header]

#Die herausfiltern
df_plot <- df_plot %>%
  filter(!(artnr %in% c(1528, 4136, 1525, 4129)))

#Zum ursprünglichen Dataframe jetzt die aggregierten hinzufügen
df_final <- rbind(as.data.frame(df_plot), as.data.frame(df_brezeln))
df_final <- rbind(as.data.frame(df_final), as.data.frame(df_laugenstange))

#keine negativen Retouren
df_final$Retoure <- if_else(df_final$endbestand < 0, 0, df_final$endbestand)

Retourenanalyse - Retouren nach Produkten & Filialen

#Auf weilerstraße filtern
pg1 <- df_final %>%
  filter(filiale == "Weilerstr.")

#Retouren aggregieren auf Artikel
pg <- aggregate(pg1$Retoure, by=list(Artikel=pg1$artnr), FUN=sum)
pg$x <- as.numeric(pg$x)
pg <- pg[order(-pg$x),] 
pg$Artikel <- as.character(pg$Artikel)
pg$Artikel <- factor(pg$Artikel, levels = c(as.character(pg$Artikel)))
colnames(pg) <- c("Artikel", "Weilerstr.")

#Dataframe zuweisen
df_pg <- pg

#Auf Rathausplatz filtern
pg1 <- df_final %>%
  filter(filiale == "Rath.Platz")
#Retouren aggregieren auf Artikel
pg <- aggregate(pg1$Retoure, by=list(Artikel=pg1$artnr), FUN=sum)
pg$x <- as.numeric(pg$x)
pg <- pg[order(-pg$x),] 
pg$Artikel <- as.character(pg$Artikel)
pg$Artikel <- factor(pg$Artikel, levels = c(as.character(pg$Artikel)))
colnames(pg) <- c("Artikel", "Rath.Platz")

#Rathausplatz Spalte hinzufügen
df_pg$Rath.Platz <- pg$Rath.Platz



plot_ly(df_pg, x = ~Artikel, y = ~Weilerstr., type = 'bar', name = 'Weilerstraße', hoverinfo = "none", marker = list(color = 'rgb(0,22,66)')) %>%
  add_trace(y = ~Rath.Platz, name = 'Rathausplatz', marker = list(color = 'rgb(158,202,225)')) %>%
  layout(xaxis = list(title = ""),
         yaxis = list(title = "", showticklabels = FALSE),
         barmode = 'group', legend = list(orientation = 'h'))

Retourenanalyse - Retourenquote nach Produkten & Filialen

library(lubridate)
df_final$jahr <- year(df_final$datum)


#Produkte und Filiale
df_retourenquote <- df_final %>% 
  group_by(artnr, filiale) %>% 
  summarise(anzahl = sum(anzahl), anfangsbestand = sum(anfangsbestand), endbestand = sum(endbestand), summe = sum(summe), std5 = sum(std5) , std6 = sum(std6), std7 = sum(std7) , std8 = sum(std8), std9 = sum(std9), std10 = sum(std10), std11 = sum(std11), std12 = sum(std12), std13 = sum(std13), std14 = sum(std14), std15 = sum(std15), std16 = sum(std16), std17 = sum(std17), std18 = sum(std18), std19 = sum(std19), std20 = sum(std20), Retoure = sum(Retoure))

#Prozentuale Retourenquote bilden
df_retourenquote$quote <- round(df_retourenquote$Retoure / df_retourenquote$anzahl,2) * 100

#Filtern nach Filialen
pg1 <- df_retourenquote %>%
  filter(filiale == "Weilerstr.")
pg2 <- df_retourenquote %>%
  filter(filiale == "Rath.Platz")

#Finalen Dataframe bilden
pg1$Rath <- pg2$quote
pg1 <- as.data.frame(pg1)
pg1$artnr <- as.character(pg1$artnr)
pg1$artnr <- factor(pg1$artnr, levels = c(as.character(pg1$artnr)))


#plot_ly(pg1, x = ~artnr, y = ~quote, type = 'bar', hoverinfo = "none", text = #paste(pg1$quote,"%",sep = ""), textposition = 'outside', name = 'Weilerstraße', marker = #list(color = 'rgb(0,22,66)')) %>%
#  add_trace(y = ~Rath, name = 'Rathausplatz', hoverinfo = "none", text = paste(pg1$Rath,"%",sep = #""), textposition = 'outside', marker = list(color = 'rgb(158,202,225)')) %>%
#  layout(xaxis = list(title = ""),
#         yaxis = list(title = "", showticklabels = FALSE),
#         barmode = 'group', legend = list(orientation = 'h'))




plot_ly(pg1, x = ~artnr, y = ~quote, type = 'bar', hoverinfo = "none", name = 'Weilerstraße', marker = list(color = 'rgb(0,22,66)')) %>%
  add_trace(y = ~Rath, name = 'Rathausplatz', hoverinfo = "none", marker = list(color = 'rgb(158,202,225)')) %>%
  layout(xaxis = list(title = ""),
         yaxis = list(title = "", showticklabels = FALSE),
         barmode = 'group', legend = list(orientation = 'h'))

Retourenanalyse - Stockoutanalyse

#Test ob Retoure oder nicht
df_final$stockout <- if_else(df_final$Retoure == 0, 1, 0)
df_final$day <- weekdays(df_final$datum)

#Initialisieren mit keinem Stockout
df_final$stockoutzeit <- 18

#Stockoutzeitpunkt festlegen -> rekursiv von letzter Öffnungszeit die Nullfolge suchen
df_final$stockoutzeit <- if_else(df_final$std17 == 0, 17, df_final$stockoutzeit)
df_final$stockoutzeit <- if_else(df_final$std17 == 0 & df_final$std16 == 0, 16, df_final$stockoutzeit)
df_final$stockoutzeit <- if_else(df_final$std17 == 0 & df_final$std16 == 0 & df_final$std15 == 0, 15, df_final$stockoutzeit)
df_final$stockoutzeit <- if_else(df_final$std17 == 0 & df_final$std16 == 0 & df_final$std15 == 0 & df_final$std14 == 0, 14, df_final$stockoutzeit)
df_final$stockoutzeit <- if_else(df_final$std17 == 0 & df_final$std16 == 0 & df_final$std15 == 0 & df_final$std14 == 0 & df_final$std13 == 0, 13, df_final$stockoutzeit)
df_final$stockoutzeit <- if_else(df_final$std17 == 0 & df_final$std16 == 0 & df_final$std15 == 0 & df_final$std14 == 0 & df_final$std13 == 0 & df_final$std12 == 0, 12, df_final$stockoutzeit)
df_final$stockoutzeit <- if_else(df_final$std17 == 0 & df_final$std16 == 0 & df_final$std15 == 0 & df_final$std14 == 0 & df_final$std13 == 0 & df_final$std12 == 0 & df_final$std11 == 0, 11, df_final$stockoutzeit)
df_final$stockoutzeit <- if_else(df_final$std17 == 0 & df_final$std16 == 0 & df_final$std15 == 0 & df_final$std14 == 0 & df_final$std13 == 0 & df_final$std12 == 0 & df_final$std11 == 0 & df_final$std10 == 0, 10, df_final$stockoutzeit)
df_final$stockoutzeit <- if_else(df_final$std17 == 0 & df_final$std16 == 0 & df_final$std15 == 0 & df_final$std14 == 0 & df_final$std13 == 0 & df_final$std12 == 0 & df_final$std11 == 0 & df_final$std10 == 0 & df_final$std9 == 0, 9, df_final$stockoutzeit)

#Nur die mit Stockout
df_stockout <- df_final %>%
  filter(stockout > 0)

#gruppieren
df_stockout <- df_stockout %>% 
  group_by(artnr, filiale) %>% 
  summarise(anzahl = sum(anzahl), anfangsbestand = sum(anfangsbestand), endbestand = sum(endbestand), summe = sum(summe), std5 = sum(std5) , std6 = sum(std6), std7 = sum(std7) , std8 = sum(std8), std9 = sum(std9), std10 = sum(std10), std11 = sum(std11), std12 = sum(std12), std13 = sum(std13), std14 = sum(std14), std15 = sum(std15), std16 = sum(std16), std17 = sum(std17), std18 = sum(std18), std19 = sum(std19), std20 = sum(std20), Retoure = sum(Retoure), stockoutzeit = mean(stockoutzeit), stockout = sum(stockout))

#Dataframe für Plot vorbereiten
pg1 <- df_stockout %>%
  filter(filiale == "Weilerstr.")
pg2 <- df_stockout %>%
  filter(filiale == "Rath.Platz")

#Durchschnittszeit
pg1$Rath <- pg2$stockout
pg1$RathZeitpunkt <- pg2$stockoutzeit
pg1$avgstockoutzeit <- (pg1$stockoutzeit + pg1$RathZeitpunkt) / 2
pg1 <- as.data.frame(pg1)


pg1$artnr <- as.character(pg1$artnr)
pg1$artnr <- factor(pg1$artnr, levels = c(as.character(pg1$artnr)))


plot_ly(pg1) %>%
  add_trace(x = ~artnr, y = ~stockout, type = 'bar', hoverinfo = "none", name = 'Weilerstraße Stockouts', marker = list(color = 'rgb(0,22,66)')) %>%
   add_trace(x = ~artnr, y = ~Rath, type = 'bar', hoverinfo = "none", name = 'Rathausplatz Stockouts', marker = list(color = 'rgb(158,202,225)')) %>%
  add_trace(x = ~artnr, y = ~avgstockoutzeit, type = 'scatter', hoverinfo = "none", mode = 'lines', name = 'Stockout Zeitpunkt', yaxis = 'y2',line = list(color = 'rgb(205, 12, 24)')) %>%
  layout(title = 'Stockoutanalyse',
         xaxis = list(title = ""),
         yaxis = list(side = 'left', title = 'Stockout Anzahl', showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
         yaxis2 = list(side = 'right', overlaying = "y", title = 'Stockout Uhrzeit'))

Data Preprocessing

Zensierte Nachfrage

Zensierte Nachfrage - Nachfragefaktoren bestimmen

#Filtern
fil <- c("Rath.Platz")
df_fil <- data %>%
  filter(filiale %in% fil)

df_fil$day <- weekdays(as.Date(df_fil$datum))

#Produktgruppe erstellen
produkt_vector <- as.data.frame(c("1525","1526","1528"))
colnames(produkt_vector) <- c("produkte")
produkt_vector$produkte <- as.character(produkt_vector$produkte)
prods <- c("1525","1526","1528")

#Wochentag Vektor erstellen
wochentag_vector <- as.data.frame(c("Montag","Dienstag","Mittwoch","Donnerstag","Freitag","Samstag","Sonntag"))
colnames(wochentag_vector) <- c("wochentage")
wochentag_vector$wochentage <- as.character(wochentag_vector$wochentage)



i <- 1
x <- 1

#Iteration über alle Produkte & Wochentage -> Nachfragefaktor getrennt nach Produkt & Wochentag
while (i <= nrow(produkt_vector)) {
  
  df_prod <- df_fil %>%
    filter(artnr == produkt_vector$produkte[i])
  x <- 1
  while(x <= nrow(wochentag_vector)) {
    
    df_sub <- df_prod %>%
      filter(day == wochentag_vector$wochentage[x])
    
    #Hier werden die Öffnungszeiten abgebildet
    if(wochentag_vector$wochentage[x] == "Samstag"){
      df_sub$isAusverkauft <- if_else(df_sub$std15 == 0, TRUE, FALSE)
    } else if(wochentag_vector$wochentage[x] == "Sonntag") {
      df_sub$isAusverkauft <- if_else(df_sub$std16 == 0, TRUE, FALSE)
    }
    else {
      df_sub$isAusverkauft <- if_else(df_sub$std17 == 0, TRUE, FALSE)
    }
    
    
    
    
    
    #1. Durchschnittswerte der Uncensored Werte bilden
    df_uc <- df_sub %>%
      filter(isAusverkauft == FALSE)
    
    #Aggregierte Nachfrage in jeder Stunde
    df_uc$H7 <- df_uc$std6 + df_uc$std7
    df_uc$H8 <- df_uc$std6 + df_uc$std7 + df_uc$std8
    df_uc$H9 <- df_uc$std6 + df_uc$std7 + df_uc$std8 + df_uc$std9
    df_uc$H10 <- df_uc$std6 + df_uc$std7 + df_uc$std8 + df_uc$std9 + df_uc$std10
    df_uc$H11 <- df_uc$std6 + df_uc$std7 + df_uc$std8 + df_uc$std9 + df_uc$std10 + df_uc$std11
    df_uc$H12 <- df_uc$std6 + df_uc$std7 + df_uc$std8 + df_uc$std9 + df_uc$std10 + df_uc$std11 + df_uc$std12
    df_uc$H13 <- df_uc$std6 + df_uc$std7 + df_uc$std8 + df_uc$std9 + df_uc$std10 + df_uc$std11 + df_uc$std12 + df_uc$std13
    df_uc$H14 <- df_uc$std6 + df_uc$std7 + df_uc$std8 + df_uc$std9 + df_uc$std10 + df_uc$std11 + df_uc$std12 + df_uc$std13 + df_uc$std14
    df_uc$H15 <- df_uc$std6 + df_uc$std7 + df_uc$std8 + df_uc$std9 + df_uc$std10 + df_uc$std11 + df_uc$std12 + df_uc$std13 + df_uc$std14 + df_uc$std15
    df_uc$H16 <- df_uc$std6 + df_uc$std7 + df_uc$std8 + df_uc$std9 + df_uc$std10 + df_uc$std11 + df_uc$std12 + df_uc$std13 + df_uc$std14 + df_uc$std15 + df_uc$std16
    df_uc$H17 <- df_uc$std6 + df_uc$std7 + df_uc$std8 + df_uc$std9 + df_uc$std10 + df_uc$std11 + df_uc$std12 + df_uc$std13 + df_uc$std14 + df_uc$std15 + df_uc$std16 + df_uc$std17
    
    df_avg_sp <- produkt_vector$produkte[i]
    df_avg_sp <- as.data.frame(df_avg_sp)
    colnames(df_avg_sp) <- c("produkt")
    
    df_avg_sp$wochentag <- wochentag_vector$wochentage[x]
    
    #Durchschnittswerte bilden
    df_avg_sp$std6 <- mean(df_uc$std6, na.rm=TRUE)
    df_avg_sp$std7 <- mean(df_uc$std7, na.rm=TRUE)
    df_avg_sp$std8 <- mean(df_uc$std8, na.rm=TRUE)
    df_avg_sp$std9 <- mean(df_uc$std9, na.rm=TRUE)
    df_avg_sp$std10 <- mean(df_uc$std10, na.rm=TRUE)
    df_avg_sp$std11 <- mean(df_uc$std11, na.rm=TRUE)
    df_avg_sp$std12 <- mean(df_uc$std12, na.rm=TRUE)
    df_avg_sp$std13 <- mean(df_uc$std13, na.rm=TRUE)
    df_avg_sp$std14 <- mean(df_uc$std14, na.rm=TRUE)
    df_avg_sp$std15 <- mean(df_uc$std15, na.rm=TRUE)
    df_avg_sp$std16 <- mean(df_uc$std16, na.rm=TRUE)
    df_avg_sp$std17 <- mean(df_uc$std17, na.rm=TRUE)
    
    df_avg_sp$H7 <- mean(df_uc$H7, na.rm=TRUE)
    df_avg_sp$H8 <- mean(df_uc$H8, na.rm=TRUE)
    df_avg_sp$H9 <- mean(df_uc$H9, na.rm=TRUE)
    df_avg_sp$H10 <- mean(df_uc$H10, na.rm=TRUE)
    df_avg_sp$H11 <- mean(df_uc$H11, na.rm=TRUE)
    df_avg_sp$H12 <- mean(df_uc$H12, na.rm=TRUE)
    df_avg_sp$H13 <- mean(df_uc$H13, na.rm=TRUE)
    df_avg_sp$H14 <- mean(df_uc$H14, na.rm=TRUE)
    df_avg_sp$H15 <- mean(df_uc$H15, na.rm=TRUE)
    df_avg_sp$H16 <- mean(df_uc$H16, na.rm=TRUE)
    df_avg_sp$H17 <- mean(df_uc$H17, na.rm=TRUE)
    
    #Rekursiv die Nachfragefaktoren bilden
    df_avg_sp$K17 <- 1
    df_avg_sp$K16 <- df_avg_sp$K17 / (1-(df_avg_sp$std17/df_avg_sp$H17))
    df_avg_sp$K15 <- df_avg_sp$K16 / (1-(df_avg_sp$std16/df_avg_sp$H16))
    df_avg_sp$K14 <- df_avg_sp$K15 / (1-(df_avg_sp$std15/df_avg_sp$H15))
    df_avg_sp$K13 <- df_avg_sp$K14 / (1-(df_avg_sp$std14/df_avg_sp$H14))
    df_avg_sp$K12 <- df_avg_sp$K13 / (1-(df_avg_sp$std13/df_avg_sp$H13))
    df_avg_sp$K11 <- df_avg_sp$K12 / (1-(df_avg_sp$std12/df_avg_sp$H12))
    df_avg_sp$K10 <- df_avg_sp$K11 / (1-(df_avg_sp$std11/df_avg_sp$H11))
    df_avg_sp$K9 <- df_avg_sp$K10 / (1-(df_avg_sp$std10/df_avg_sp$H10))
    df_avg_sp$K8 <- df_avg_sp$K9 / (1-(df_avg_sp$std9/df_avg_sp$H9))
    df_avg_sp$K7 <- df_avg_sp$K8 / (1-(df_avg_sp$std8/df_avg_sp$H8))
    df_avg_sp$K6 <- df_avg_sp$K7 / (1-(df_avg_sp$std7/df_avg_sp$H7))
    
    if(x == 1 & i == 1) {
      df_avg <- df_avg_sp  
    } else {
      df_avg <- rbind(df_avg, df_avg_sp)
    }
    
    x <- x +1
  }
  
  i <- i +1  
}

Zensierte Nachfrage - Stockoutzeitpunkt bestimmen

#Spalten selektieren und runden
df_arbeit <- df_avg %>%
  select(produkt, wochentag, K6,K7,K8,K9,K10,K11,K12,K13,K14,K15,K16,K17)
df_arbeit$K17 <- round(df_arbeit$K17,2)
df_arbeit$K16 <- round(df_arbeit$K16,2)
df_arbeit$K15 <- round(df_arbeit$K15,2)
df_arbeit$K14 <- round(df_arbeit$K14,2)
df_arbeit$K13 <- round(df_arbeit$K13,2)
df_arbeit$K12 <- round(df_arbeit$K12,2)
df_arbeit$K11 <- round(df_arbeit$K11,2)
df_arbeit$K10 <- round(df_arbeit$K10,2)
df_arbeit$K9 <- round(df_arbeit$K9,2)
df_arbeit$K8 <- round(df_arbeit$K8,2)
df_arbeit$K7 <- round(df_arbeit$K7,2)
df_arbeit$K6 <- round(df_arbeit$K6,2)


#Stockoutzeitpunkt berechnen
df_prod <- df_fil %>%
  filter(artnr %in% prods)

#Öffnungszeiten abbilden
df_prod$isAusverkauft <- if_else(df_prod$std17 == 0 & df_prod$day != "Samstag" & df_prod$day != "Sonntag", TRUE, FALSE)
df_prod$isAusverkauft <- if_else(df_prod$std15 == 0 & df_prod$day == "Samstag", TRUE, df_prod$isAusverkauft)
df_prod$isAusverkauft <- if_else(df_prod$std16 == 0 & df_prod$day == "Sonntag", TRUE, df_prod$isAusverkauft)

#Wenn Ausverkauft ist -> zensiert
df_c_sub <- df_prod %>%
  filter(isAusverkauft == TRUE)

#Rekursiv werden Zeitpunkte der zensierten Nachfrage bestimmt
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0, 16, 17)
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0 & df_c_sub$std15 == 0, 15, df_c_sub$StockoutZeitpunkt)
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0 & df_c_sub$std15 == 0 & df_c_sub$std14 == 0, 14, df_c_sub$StockoutZeitpunkt)
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0 & df_c_sub$std15 == 0 & df_c_sub$std14 == 0 & df_c_sub$std13 == 0, 13, df_c_sub$StockoutZeitpunkt)
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0 & df_c_sub$std15 == 0 & df_c_sub$std14 == 0 & df_c_sub$std13 == 0 & df_c_sub$std12 == 0, 12, df_c_sub$StockoutZeitpunkt)
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0 & df_c_sub$std15 == 0 & df_c_sub$std14 == 0 & df_c_sub$std13 == 0 & df_c_sub$std12 == 0 & df_c_sub$std11 == 0, 11, df_c_sub$StockoutZeitpunkt)
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0 & df_c_sub$std15 == 0 & df_c_sub$std14 == 0 & df_c_sub$std13 == 0 & df_c_sub$std12 == 0 & df_c_sub$std11 == 0 & df_c_sub$std10 == 0, 10, df_c_sub$StockoutZeitpunkt)
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0 & df_c_sub$std15 == 0 & df_c_sub$std14 == 0 & df_c_sub$std13 == 0 & df_c_sub$std12 == 0 & df_c_sub$std11 == 0 & df_c_sub$std10 == 0 & df_c_sub$std9 == 0, 9, df_c_sub$StockoutZeitpunkt)
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0 & df_c_sub$std15 == 0 & df_c_sub$std14 == 0 & df_c_sub$std13 == 0 & df_c_sub$std12 == 0 & df_c_sub$std11 == 0 & df_c_sub$std10 == 0 & df_c_sub$std9 == 0 & df_c_sub$std8 == 0, 8, df_c_sub$StockoutZeitpunkt)
df_c_sub$StockoutZeitpunkt <- if_else(df_c_sub$std16 == 0 & df_c_sub$std15 == 0 & df_c_sub$std14 == 0 & df_c_sub$std13 == 0 & df_c_sub$std12 == 0 & df_c_sub$std11 == 0 & df_c_sub$std10 == 0 & df_c_sub$std9 == 0 & df_c_sub$std8 == 0 & df_c_sub$std7 == 0, 7, df_c_sub$StockoutZeitpunkt)

test <- df_c_sub %>% 
  group_by(datum) %>% 
  summarise(MinStockoutZeitpunkt = min(StockoutZeitpunkt))

df_final <- left_join(df_c_sub, test, by = "datum")

Zensierte Nachfrage - Reparieren der zensierten Nachfrage

#Dummy Werte
df_prod$StockoutZeitpunkt <- 500
df_prod$MinStockoutZeitpunkt <- 500

#Alle unter die mit Stockout binden
df_final <- rbind(df_final, df_prod)

#Die doppelten raus -> es wird immer das nächste duplikat entfernt
df_final  <- df_final[!duplicated(df_final[2:4]),]

#Hier wird jetzt auch für "vollständige" Werte die keinen Stockout erfahren haben der Minimale Stockoutzeitpunkt hinzugefügt für die Berechnung der Substitutionseffekte
test <- df_final %>%
  group_by(datum) %>%
  summarise(MinStockoutZeitpunktNeu = min(StockoutZeitpunkt))

df_final <- left_join(df_final, test, by = "datum")

#Die Werte enthalten den kompletten Demand und sind somit vollständig
df_complete_uncensored <- df_final %>%
  filter(MinStockoutZeitpunktNeu == 500)

#Diese werden rausgefiltert
df_final <- df_final %>%
  filter(MinStockoutZeitpunktNeu < 500)


#Dummy Werte
df_final$demand <- 0
df_final$censoredDemand <- 0

#Iteration über jede Zeile
i <- 1
while (i <= nrow(df_final)) {
  
  #Ich hol mir den Produktname und den Tagname damit ich in meiner avg tabelle den richtigen wert bekomme  
  produkt_name <- df_final$artnr[i]
  Tag_name <- df_final$day[i]
  
  df_avg_loop <- df_avg %>%
    filter(produkt == produkt_name) %>%
    filter(wochentag == Tag_name)
  
  #Hier tatsächlich den echten Wert anhand des MinStockoutZeitpunktes
  #Je nach Zeitpunkt wird der Demand nur bis zu diesem Zeitpunkt genommen
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 17, ((df_avg_loop$K17 + df_avg_loop$K16)/2)*rowSums(df_final[, c(10:21)])[i],0),0)
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 16, ((df_avg_loop$K16 + df_avg_loop$K15)/2)*rowSums(df_final[, c(10:20)])[i],df_final$demand[i]),0)
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 15, ((df_avg_loop$K15 + df_avg_loop$K14)/2)*rowSums(df_final[, c(10:19)])[i],df_final$demand[i]),0)
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 14, ((df_avg_loop$K14 + df_avg_loop$K13)/2)*rowSums(df_final[, c(10:18)])[i],df_final$demand[i]),0)
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 13, ((df_avg_loop$K13 + df_avg_loop$K12)/2)*rowSums(df_final[, c(10:17)])[i],df_final$demand[i]),0)
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 12, ((df_avg_loop$K12 + df_avg_loop$K11)/2)*rowSums(df_final[, c(10:16)])[i],df_final$demand[i]),0)
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 11, ((df_avg_loop$K11 + df_avg_loop$K10)/2)*rowSums(df_final[, c(10:15)])[i],df_final$demand[i]),0)
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 10, ((df_avg_loop$K10 + df_avg_loop$K9)/2)*rowSums(df_final[, c(10:14)])[i],df_final$demand[i]),0)
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 9, ((df_avg_loop$K9 + df_avg_loop$K8)/2)*rowSums(df_final[, c(10:13)])[i],df_final$demand[i]),0)
  df_final$demand[i] <- round(if_else(df_final$MinStockoutZeitpunktNeu[i] == 8, ((df_avg_loop$K8 + df_avg_loop$K7)/2)*rowSums(df_final[, c(10:12)])[i],df_final$demand[i]),0)
  
  #Hier nur den Wert wenn ich rein die Censored Daten ohne Substitution betrachte
  #Ich nehme also nicht den Minstockoutzeitpunkt sondern nur den Stockoutzeitpunkt
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 17, ((df_avg_loop$K17 + df_avg_loop$K16)/2)*rowSums(df_final[, c(10:21)])[i],0),0)
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 16, ((df_avg_loop$K16 + df_avg_loop$K15)/2)*rowSums(df_final[, c(10:20)])[i],df_final$censoredDemand[i]),0)
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 15, ((df_avg_loop$K15 + df_avg_loop$K14)/2)*rowSums(df_final[, c(10:19)])[i],df_final$censoredDemand[i]),0)
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 14, ((df_avg_loop$K14 + df_avg_loop$K13)/2)*rowSums(df_final[, c(10:18)])[i],df_final$censoredDemand[i]),0)
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 13, ((df_avg_loop$K13 + df_avg_loop$K12)/2)*rowSums(df_final[, c(10:17)])[i],df_final$censoredDemand[i]),0)
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 12, ((df_avg_loop$K12 + df_avg_loop$K11)/2)*rowSums(df_final[, c(10:16)])[i],df_final$censoredDemand[i]),0)
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 11, ((df_avg_loop$K11 + df_avg_loop$K10)/2)*rowSums(df_final[, c(10:15)])[i],df_final$censoredDemand[i]),0)
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 10, ((df_avg_loop$K10 + df_avg_loop$K9)/2)*rowSums(df_final[, c(10:14)])[i],df_final$censoredDemand[i]),0)
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 9, ((df_avg_loop$K9 + df_avg_loop$K8)/2)*rowSums(df_final[, c(10:13)])[i],df_final$censoredDemand[i]),0)
  df_final$censoredDemand[i] <- round(if_else(df_final$StockoutZeitpunkt[i] == 8, ((df_avg_loop$K8 + df_avg_loop$K7)/2)*rowSums(df_final[, c(10:12)])[i],df_final$censoredDemand[i]),0)
  
  i <- i +1
}

Zensierte Nachfrage - Fallbetrachtung bei 3 Produkten

Hier muss überprüft werden wie hoch die Substitution ist & wie die Substitution verteilt werden kann auf die Produkte.

#zeitlich ordnen
df_final <- df_final[order(as.Date(df_final$datum)),]

#Wenn ausverkauft ist, censored Demand nehmen und ansonsten Actuals
df_final$censoredDemand <- if_else(df_final$isAusverkauft == FALSE, df_final$anzahl, df_final$censoredDemand)


#Herausfinden von des Produktes mit MinimalZeitpunkt -> Ausgangslage für Substitution
df_final$CensoredDemandZuwachs <- 0

#Veränderung von Demand Prognose vs. Actuals
df_final$CensoredDemandZuwachs <- df_final$demand - df_final$anzahl

#Wenn Veränderung > 0 bei isAusverkauft == FALSE -> dann auf 0 setzen -> Modell sagt mehr als tatsächlich verkauft -> geht nicht
df_final$CensoredDemandZuwachs <- if_else(df_final$isAusverkauft == FALSE & df_final$CensoredDemandZuwachs > 0, 0, df_final$CensoredDemandZuwachs)


#Da wir mehrere Produkte haben muss man die Fallunterscheidungen immer an einem Tag in der bestimmten Gruppe betrachten
#Deswegen werden zunächst die Daten gruppiert nach isAusverkauft und datum -> pro Datum 2 Werte mit ausverkauft und nicht ausverkauft

#Jetzt die Werte gruppieren, um die Gesamtveränderung zu bekommen
test <- df_final %>%
  group_by(datum, isAusverkauft) %>%
  summarise(GesamtVeraenderung = sum(CensoredDemandZuwachs))

df_final <- left_join(df_final, test, by = c("datum" = "datum", "isAusverkauft" = "isAusverkauft"))

#Die Spalte Gesamtveränderung zeigt für jedes Datum und jede isAusverkauft Variante an wie viel insgesamt -> anhand dessen mach ich die Fallunterscheidung


#Fall1: Idealfall

#Gesamtveränderung -> hier schau ich mir jetzt insgesamt an ob die Gesamtveränderung > 0 ist.
Case <- test %>%
  group_by(datum) %>%
  summarise(case = sum(GesamtVeraenderung))

#Wenn True, dann bewegt sich der Bereich zwischen 0 und der Veränderung -> Also wenn TRUE sind wir im Fall1 (Idealfall)
Case$CaseRichtig <- if_else(Case$case > 0, TRUE, FALSE)

#wieder an test joinen -> jetzt weiß ich, dass CaseRichtig == TRUE mein Idealfall ist
test <- left_join(test, Case, by = "datum")

#Da wo CaseRichtig und isAusverkauft == TRUE handelt es sich um den Idealfall -> danach filter ich
df_richtig <- test %>%
  filter(isAusverkauft == TRUE) %>%
  filter(CaseRichtig == TRUE) %>%
  select(datum, GesamtVeraenderung)
colnames(df_richtig) <- c("datum","Teiler")

#In der Spalte Teiler steht jetzt immer der Wert an zusätzlichem Censored Demand von dem ersten Produkt
#Eigentlich brauche ich das nicht wirklich weil der Predictete Demand ja bereits korrekt ist. Ich mach es trotzdem um später die Substiutionsrate zu berechnen

#Diese Daten dran joinen -> jetzt kann im Idealfall einfach die Substitutionsrate berechnet werden
df_final <- left_join(df_final, df_richtig, by = "datum")




#Fall2: Modell sagt zu wenig -> 100% Substitution

#Es kann maximal die Menge reduziert werden, die das Stockoutprodukt hinzubekommt (Upper Bound)
#Deswegen wird die Gesamtmenge der Reduktion genommen und die Reduktionsmengen der einzelnen Produkte damit geteilt -> Dadurch erhält man prozentuale Werte
#Mit diesem prozentualen Wert kann die Upper Bound Menge multipliziert werden, sodass sie gerecht aufgeteilt wird.
#Für Fall 2 = Modell sagt zu wenig -> Ich filter hier einfach nach CaseRichtig == FALSE -> Also nicht Idealfall
df_fall2 <- test %>%
  filter(isAusverkauft == TRUE) %>%
  filter(CaseRichtig == FALSE) %>%
  select(datum, GesamtVeraenderung)
colnames(df_fall2) <- c("datum","TeilerFall2")

#Hier gibt es jetzt nochmal einen Fallunterschied. Dies ist Fall3. Bei Fall3 sind mehrere Produkte ausverkauft, und diese sind in Summe nicht Fall 1 zuzuordnen.
#Dieser Fall muss wieder extra betrachtet werden

#Brauch ich für später unten für Fall 3 -> Modell sagt insgesamt zu wenig
df_fall3 <- df_fall2 %>%
  filter(TeilerFall2 < 0)

#Die restlichen rausfiltern -> Fall 2
df_fall2 <- df_fall2 %>%
  filter(TeilerFall2 >= 0)

df_final <- left_join(df_final, df_fall2, by = "datum")
#In TeilerFall2 sehe ich jetzt den zusätzlichen Censored Demand. Diesen benötige ich, um meine Substiutionsrate mit dem Teiler zu multiplizieren

Zensierte Nachfrage - Substitutionsfaktoren

#Substitutionsrate für Fall 1 -> brauche ich eigentlich nicht zu berechnen, um den Demand zu bekommen
df_final$SubstitutsrateFall1 <- 0
df_final$SubstitutsrateFall1 <- if_else(df_final$isAusverkauft == FALSE, round((df_final$CensoredDemandZuwachs / df_final$Teiler)*(-1),2), 0)

#Substituonsrate für Fall2 -> die brauche ich. In der Spalte Gesamtveränderung sehe ich die summierten Werte für isAusverkauft == FALSE. 
#Also nehme ich Die Gesamtmenge die reduziert werden soll und Teile dann die Menge die reduziert werden soll von Produkt i mit der Gesamtmenge
df_final$SubstitutsrateFall2 <- 0
df_final$SubstitutsrateFall2 <- if_else(df_final$isAusverkauft == FALSE, round((df_final$CensoredDemandZuwachs / df_final$GesamtVeraenderung),2), 0)

#Das ist die Menge, die abgezogen werden muss im Fall 2. 
df_final$Fall2DemandVerbesserung <- df_final$TeilerFall2 * df_final$SubstitutsrateFall2

#Sonst überall 0 abziehen
df_final$Fall2DemandVerbesserung[is.na(df_final$Fall2DemandVerbesserung)] <- 0

#Wenn isAusverkauft == FALSE die Menge abziehen, um den neuen "echten" Demand zu bekommen
df_final$FinalDemandFall2 <- if_else(df_final$isAusverkauft == FALSE, df_final$anzahl - df_final$Fall2DemandVerbesserung, df_final$demand)
df_final$FinalDemandFall2 <- round(df_final$FinalDemandFall2,0) 


#Fall3: Mehrere Produkte haben Stockout und Modell sagt insgesamt zu wenig
#Trifft weder auf Fall1 noch Fall2 zu -> Seperat betrachten

#Ich hol mir die entsprechenden Daten
fall3_vec <- df_fall3$datum

#Das sind jetzt alle Daten, bei denen Fall 3 der Fall ist: Fall 3 ist wenn mehrere Stockouts sind und das Modell insgesamt wiederzu wenig predictet
df_fall3 <- df_final %>%
  filter(datum %in% fall3_vec)

#Das Max raussuchen -> Das Max ist die censored Menge, die hinzukommt. Da beim anderen Produkt (das out of Stock ist), die predictete Menge unterhalb der Actuals ist, gibt es keinen zusätzlicehn Demand
#Diese Menge muss jetzt wieder verteilt werden -> Logik wie in Fall2
fall3 <- df_fall3 %>%
  group_by(datum) %>%
  summarise(Max = max(CensoredDemandZuwachs))

#Herausfiltern von negativen -> Das kommt wenn bspw. Verkäuferinnen nach Ladenschluss noch was holen und das wird dann erst unter 18 Uhr gespiechert
#Sollte eigentlich nicht vorkommen, weil oben die Werte genullt werden
fall3 <- fall3 %>%
  filter(Max > 0)

#Herausfinden, wo der Demandzuwachs < 0 ist -> Das heißt Modell hat zu wenig predictet
df_fall3$isNegative <- if_else(df_fall3$CensoredDemandZuwachs < 0, TRUE, FALSE)

#Diese negativen werden summiert als Basis um den Faktor zu berechnen
df_fall3Teiler <- df_fall3 %>%
  group_by(datum, isNegative) %>%
  summarise(Teiler = sum(CensoredDemandZuwachs))

df_fall3Teiler <- df_fall3Teiler %>%
  filter(Teiler < 0)
#Wieder die Nuller rausfiltern

#Jetzt die Zusammenjoinen -> Also hab ich in Fall3Max die Menge, die verteilt werden muss und in Fall3Teiler die Menge mit der ich den Faktor berechnen kann
fall3 <- left_join(fall3, df_fall3Teiler, by = "datum")
fall3 <- fall3 %>%
  select(-isNegative)
colnames(fall3) <- c("datum","Fall3Max","Fall3Teiler")

#An den finalen ran joinen
df_final <- left_join(df_final, fall3, by = "datum")
#Das soll nur gemacht werden, für die Werte mit Substitution -> Menge die von den Actuals abgezogen werden muss.
df_final$Fall3DemandVerbesserung <- if_else(df_final$CensoredDemandZuwachs != df_final$Fall3Max, (df_final$CensoredDemandZuwachs / df_final$Fall3Teiler) * df_final$Fall3Max, 0)
#Rest sind Nuller
df_final$Fall3DemandVerbesserung[is.na(df_final$Fall3DemandVerbesserung)] <- 0

#Substitutionsrate Fall3
df_final$SubstitutsrateFall3 <- if_else(df_final$CensoredDemandZuwachs != df_final$Fall3Max, round(df_final$CensoredDemandZuwachs / df_final$Fall3Teiler,2), 0)




#Finaler Demand über alle Fälle berechnen
#Fall1: demand ist richtig
#Fall2: Actuals - Fall2Demandverbesserung
#Fall3: Actuals - Fall3Demandverbesserung

df_final$FinalerDemand <- df_final$demand

#Fall2:
df_final$FinalerDemand <- if_else(is.na(df_final$TeilerFall2) == FALSE & df_final$TeilerFall2 != df_final$CensoredDemandZuwachs, df_final$anzahl - df_final$Fall2DemandVerbesserung, df_final$FinalerDemand)

#Fall3:
df_final$FinalerDemand <- if_else(is.na(df_final$Fall3Teiler) == FALSE & df_final$Fall3Max != df_final$CensoredDemandZuwachs, df_final$anzahl - df_final$Fall3DemandVerbesserung, df_final$FinalerDemand)

#Runden
df_final$FinalerDemand <- round(df_final$FinalerDemand,0)



#Finale Substitutionsraten berechnen
df_final$FinaleSubstitutionsrate <- 0

#Fall1
df_final$FinaleSubstitutionsrate <- if_else(is.na(df_final$Teiler) == FALSE & df_final$Teiler != df_final$CensoredDemandZuwachs, df_final$SubstitutsrateFall1, df_final$FinaleSubstitutionsrate)
df_final$FinaleSubstitutionsrate <- if_else(is.na(df_final$Teiler) == FALSE & df_final$Teiler == df_final$GesamtVeraenderung, 500, df_final$FinaleSubstitutionsrate)

#die 500 sind ein Dummywert damit ich das später rausfiltern kann -> Das ist ja immer das ausgehende Produkt

#Fall2
df_final$FinaleSubstitutionsrate <- if_else(is.na(df_final$TeilerFall2) == FALSE & df_final$TeilerFall2 != df_final$CensoredDemandZuwachs, df_final$SubstitutsrateFall2, df_final$FinaleSubstitutionsrate)
df_final$FinaleSubstitutionsrate <- if_else(is.na(df_final$TeilerFall2) == FALSE & df_final$TeilerFall2 == df_final$GesamtVeraenderung, 500, df_final$FinaleSubstitutionsrate)

#Fall3
df_final$FinaleSubstitutionsrate <- if_else(is.na(df_final$Fall3Teiler) == FALSE & df_final$Fall3Max != df_final$CensoredDemandZuwachs, df_final$SubstitutsrateFall3, df_final$FinaleSubstitutionsrate)
df_final$FinaleSubstitutionsrate <- if_else(is.na(df_final$Fall3Teiler) == FALSE & df_final$Fall3Max == df_final$CensoredDemandZuwachs, 500, df_final$FinaleSubstitutionsrate)


#Substitutionsrate berechnen für Brezeln
#Nur wenn die Brezel das einzig ausverkaufte Produkt war, kann die Substitutionsrate berechnet werden
df_brezel <- df_final %>%
  filter(artnr == "1525") %>%
  filter(isAusverkauft == TRUE) %>%
  filter(CensoredDemandZuwachs == GesamtVeraenderung)
brezel_vec <- df_brezel$datum

df_brezel <- df_final %>%
  filter(datum %in% brezel_vec) %>%
  filter(FinaleSubstitutionsrate < 500)

brezel_laugenstange <- df_brezel %>%
  filter(artnr == "1528")
Sub_brezel_stange <- mean(brezel_laugenstange$FinaleSubstitutionsrate)

brezel_laugenbroetchen <- df_brezel %>%
  filter(artnr == "1526")
Sub_brezel_broetchen <- mean(brezel_laugenbroetchen$FinaleSubstitutionsrate)

Sub_brezel_broetchen + Sub_brezel_stange

## [1] 0.7453111

#Substitutionsrate berechnen für Laugenstangen
df_brezel <- df_final %>%
  filter(artnr == "1528") %>%
  filter(isAusverkauft == TRUE) %>%
  filter(CensoredDemandZuwachs == GesamtVeraenderung)
brezel_vec <- df_brezel$datum

df_brezel <- df_final %>%
  filter(datum %in% brezel_vec) %>%
  filter(FinaleSubstitutionsrate < 500)

laugenstange_brezel <- df_brezel %>%
  filter(artnr == "1525")
Sub_stange_brezel <- mean(laugenstange_brezel$FinaleSubstitutionsrate)

laugenstange_laugenbroetchen <- df_brezel %>%
  filter(artnr == "1526")
Sub_stange_broetchen <- mean(laugenstange_laugenbroetchen$FinaleSubstitutionsrate)

Sub_stange_broetchen + Sub_stange_brezel

## [1] 0.6739878

#Substitutionsrate berechnen für Laugenbrötchen
df_brezel <- df_final %>%
  filter(artnr == "1526") %>%
  filter(isAusverkauft == TRUE) %>%
  filter(CensoredDemandZuwachs == GesamtVeraenderung)
brezel_vec <- df_brezel$datum

df_brezel <- df_final %>%
  filter(datum %in% brezel_vec) %>%
  filter(FinaleSubstitutionsrate < 500)

broetchen_brezel <- df_brezel %>%
  filter(artnr == "1525")
Sub_broetchene_brezel <- mean(broetchen_brezel$FinaleSubstitutionsrate)

laugenbroetchen_laugenstange <- df_brezel %>%
  filter(artnr == "1528")
Sub_broetchene_stange <- mean(laugenbroetchen_laugenstange$FinaleSubstitutionsrate)

Sub_broetchene_brezel + Sub_broetchene_stange

## [1] 0.6120143

round(Sub_brezel_broetchen,4)

## [1] 0.3792

round(Sub_brezel_stange,4)

## [1] 0.3661

round(Sub_stange_brezel,4)

## [1] 0.4035

round(Sub_stange_broetchen,4)

## [1] 0.2705

round(Sub_broetchene_brezel,4)

## [1] 0.3102

round(Sub_broetchene_stange,4)

## [1] 0.3018

df <- df_final %>%
  select(datum, filiale, artnr, anfangsbestand, FinalerDemand, X)
colnames(df) <- c("Datum", "KundenNr", "ArtNr", "GelieferteMenge", "VerkaufteMenge", "Artikelgruppe")

Zensierte Nachfrage - Darstellung der Substitutionsfaktoren

df_test <- df_final %>%
  filter(day == "Dienstag")

df_p1 <- df_test %>%
  filter(artnr == 1525)
df_p2 <- df_test %>%
  filter(artnr == 1528)
df_p3 <- df_test %>%
  filter(artnr == 1526)





ax <- list(
  showticklabels = FALSE,
  showgrid = FALSE
)


p1 <- plot_ly(df_p1, x = ~datum, y = ~FinalerDemand, hoverinfo = "none", name = 'demand', type = 'scatter', mode = 'lines', line = list(color = 'transparent'),showlegend = FALSE) %>%
  add_trace(y = ~anzahl, hoverinfo = "none", name = 'Lower', type = 'scatter', mode = 'lines', line = list(color = 'transparent'),showlegend = FALSE, fill = 'tonexty', fillcolor='rgba(31, 119, 180,0.2)') %>%
  add_trace(y = ~FinalerDemand, hoverinfo = "none", name = 'demand_repariert', type = 'scatter', mode = 'lines', line = list(color = 'rgb(8,48,107)')) %>%
  add_trace(y = ~anzahl,hoverinfo = "none", name = 'Actuals', type = 'scatter', mode = 'lines', line = list(color = 'rgb(205, 12, 24)')) %>%
  layout(xaxis = ax, yaxis = list(showticklabels = FALSE))
p2 <- plot_ly(df_p2, x = ~datum, y = ~FinalerDemand, hoverinfo = "none", name = 'demand', type = 'scatter', mode = 'lines', line = list(color = 'transparent'),showlegend = FALSE) %>%
  add_trace(y = ~anzahl, hoverinfo = "none", name = 'Lower', type = 'scatter', mode = 'lines', line = list(color = 'transparent'),showlegend = FALSE, fill = 'tonexty', fillcolor='rgba(31, 119, 180,0.2)') %>%
  add_trace(y = ~FinalerDemand, hoverinfo = "none", name = 'demand_repariert', type = 'scatter', mode = 'lines', line = list(color = 'rgb(8,48,107)')) %>%
  add_trace(y = ~anzahl, hoverinfo = "none", name = 'Actuals', type = 'scatter', mode = 'lines', line = list(color = 'rgb(205, 12, 24)')) %>%
  layout(xaxis = ax, yaxis = list(showticklabels = FALSE))
p3 <- plot_ly(df_p3, x = ~datum, y = ~FinalerDemand, hoverinfo = "none", name = 'demand', type = 'scatter', mode = 'lines', line = list(color = 'transparent'),showlegend = FALSE) %>%
  add_trace(y = ~anzahl, hoverinfo = "none", name = 'Lower', type = 'scatter', mode = 'lines', line = list(color = 'transparent'),showlegend = FALSE, fill = 'tonexty', fillcolor='rgba(31, 119, 180,0.2)') %>%
  add_trace(y = ~FinalerDemand, hoverinfo = "none", name = 'demand_repariert', type = 'scatter', mode = 'lines', line = list(color = 'rgb(8,48,107)')) %>%
  add_trace(y = ~anzahl, hoverinfo = "none", name = 'Actuals', type = 'scatter', mode = 'lines', line = list(color = 'rgb(205, 12, 24)'))%>%
  layout(xaxis = ax, yaxis = list(showticklabels = FALSE))
p <- subplot(p1, p2, p3, nrows = 3)
p

Feature Engineering

Ferien Features

library(httr)
library(jsonlite)

#------------------- 1. Ferien Termine-----------------

ferien <- GET("https://ferien-api.de/api/v1/holidays/BW")
#API Callen
ferien_text <- content(ferien, "text")
#Text bekommen
ferien_json <- fromJSON(ferien_text, flatten = FALSE)
#Json File lesen
ferien_df <- as.data.frame(ferien_json)
#Als Dataframe speichern

ferien_df$startdatum <- as.Date(as.POSIXct(ferien_df$start ,format = "%Y-%m-%dT%H:%M"))
ferien_df$enddatum <- as.Date(as.POSIXct(ferien_df$end ,format = "%Y-%m-%dT%H:%M"))
#Datum speichern

#Für jede Ferien wird die jeweilige Spalte übergeben, sodass ich einen Dataframe bekomme mit dem jeweiligen Ferientag für jedes Datum
holidays <- data.frame(matrix(ncol = 4, nrow = 0))
colnames(holidays) <- c("Holiday","Datum","Lower","Upper")

i <- 1
while(i <= nrow(ferien_df)) {
  df_hol <- data_frame(
    Holiday = ferien_df[i,5],
    Datum = seq(as.Date(ferien_df[i,7]), as.Date(ferien_df[i,8]), "days"),
    Lower = -1,
    Upper = 1
  )
  
  holidays <- rbind(holidays, df_hol)
  
  i <- i + 1
}




#2. -------------------- Feiertage Termine--------------------

#2017
ferien <- GET("https://feiertage-api.de/api/?jahr=2017&nur_land=BW")
#Get Ferien Termine
ferien_text <- content(ferien, "text")
ferien_json <- fromJSON(ferien_text, flatten = FALSE)
#Json als List speichern
df_feiertage <- data.frame(matrix(unlist(ferien_json), nrow=14, byrow=T),stringsAsFactors=FALSE)
#Speichert die List inhalte in einen Dataframe
df_feiertage$feiertag_name <- names(ferien_json)
#Die Namen der ersten Liste sind die Feiertagsnamen

df_feiertage <- df_feiertage %>%
  select(feiertag_name, X1)
colnames(df_feiertage) <- c("Holiday", "Datum")
df_feiertage$Lower <- -1
df_feiertage$Upper <- 1

holidays <- rbind(holidays, df_feiertage)


#2018
ferien <- GET("https://feiertage-api.de/api/?jahr=2018&nur_land=BW")
#Get Ferien Termine
ferien_text <- content(ferien, "text")
ferien_json <- fromJSON(ferien_text, flatten = FALSE)
#Json als List speichern
df_feiertage <- data.frame(matrix(unlist(ferien_json), nrow=14, byrow=T),stringsAsFactors=FALSE)
#Speichert die List inhalte in einen Dataframe
df_feiertage$feiertag_name <- names(ferien_json)
#Die Namen der ersten Liste sind die Feiertagsnamen

df_feiertage <- df_feiertage %>%
  select(feiertag_name, X1)
colnames(df_feiertage) <- c("Holiday", "Datum")
df_feiertage$Lower <- -1
df_feiertage$Upper <- 1

holidays <- rbind(holidays, df_feiertage)


#2019
ferien <- GET("https://feiertage-api.de/api/?jahr=2019&nur_land=BW")
#Get Ferien Termine
ferien_text <- content(ferien, "text")
ferien_json <- fromJSON(ferien_text, flatten = FALSE)
#Json als List speichern
df_feiertage <- data.frame(matrix(unlist(ferien_json), nrow=14, byrow=T),stringsAsFactors=FALSE)
#Speichert die List inhalte in einen Dataframe
df_feiertage$feiertag_name <- names(ferien_json)
#Die Namen der ersten Liste sind die Feiertagsnamen

df_feiertage <- df_feiertage %>%
  select(feiertag_name, X1)
colnames(df_feiertage) <- c("Holiday", "Datum")
df_feiertage$Lower <- -1
df_feiertage$Upper <- 1

holidays <- rbind(holidays, df_feiertage)
colnames(holidays) <- c("holiday","ds","lower_window","upper_window")

holidays_ml <- holidays %>%
  select(ds, holiday)
colnames(holidays_ml) <- c("Datum","holiday")

df <- left_join(df, holidays_ml, by = "Datum")
remove(holidays_ml)

Wetter Feature

library("rdwd")

stations <- nearbyStations(48.165531, 8.6251283, 20)
#Hier sind die verfügbaren Stationen mit den entsprechend verfügbaren Variablen

#Daten werden gesplittet in Recent und Historische -> Beide abrufen und anschließend joinen

#Gleiches für andere Datenart
link <- selectDWD("Rottweil", res="daily", var="kl", per="recent")
file <- dataDWD(link, read=FALSE, dir="DWDdata", quiet=TRUE, force=NA)
kl <- readDWD(file, varnames=TRUE)


#Gleiches für andere Datenart
link <- selectDWD("Rottweil", res="daily", var="kl", per="historical")
file <- dataDWD(link, read=FALSE, dir="DWDdata", quiet=TRUE, force=NA)
kl_hist <- readDWD(file, varnames=TRUE)


#-------------------------------2. Forecast Wetter Daten mit Opean Weather
#Gibt für die nächsten 5 Tage die Wetterdaten aggregiert auf 3 Stunden Basis
library(owmr)
citation("owmr")


owmr_settings("ccf39b0a06683bfee02715cef67146a9")
OWM_API_KEY <- "ccf39b0a06683bfee02715cef67146a9"
Sys.setenv(OWM_API_KEY = "ccf39b0a06683bfee02715cef67146a9")
#API Keys setzen

forecast <- get_forecast("Rottweil,GER", units = "metric") %>%
  owmr_as_tibble()

library(lubridate)
library(tidyverse)

forecast$dt_txt <- as.POSIXct(forecast$dt_txt ,format = "%Y-%m-%d  %H:%M:%OS")
#Als Postixct speichern
forecast$datum <- as.Date(forecast$dt_txt)
#Uhrzeit weg

forecast_weather <- forecast %>%
  group_by(datum) %>%
  summarise(temp = mean(temp, na.rm = TRUE), temp_max = max(temp_max, na.rm = TRUE), temp_min = min(temp_min, na.rm = TRUE), feuchtigkeit = mean(humidity, na.rm = TRUE), niederschlag = sum(rain_3h, na.rm = TRUE))
#Auf Tagesbasis aggregieren

colnames(forecast_weather) <- c("Datum","Temperatur","Max_Temperatur","Min_Temperatur","Luftfeuchtigkeit","Niederschlag")


#----------------3. historische und zukünftige Werte joinen ------------

#---Historische und aktuelle joinen
df_weather <- rbind(kl_hist, kl)
df_weather <- unique(df_weather)
#Aus API des deutschen Wetterdienstes historische mit recent Daten anbinden und die doppelten Werte raus

df_weather <- df_weather %>%
  select(MESS_DATUM, TMK.Lufttemperatur, TXK.Lufttemperatur_Max, TNK.Lufttemperatur_Min, UPM.Relative_Feuchte, RSK.Niederschlagshoehe)
colnames(df_weather) <- c("Datum","Temperatur","Max_Temperatur","Min_Temperatur","Luftfeuchtigkeit","Niederschlag")

df_weather <- rbind(df_weather, forecast_weather)
#Forecast und vergangene Daten joinen


#Grundidee: Verkäufe hängen nicht vom genauen Wetter ab, sondern von Wettergruppen
#Ist es besonders warm oder besonders kalt

library(Hmisc)
df_weather <- na.omit(df_weather)

#1. Temperatur
plot(density(df_weather$Temperatur))
summary(df_weather$Temperatur)

df_weather$Temperatur_Groups <- as.numeric(cut2(df_weather$Temperatur, g=8))
df_weather$MaxTemperatur_Groups <- as.numeric(cut2(df_weather$Max_Temperatur, g=8))
df_weather$MinTemperatur_Groups <- as.numeric(cut2(df_weather$Min_Temperatur, g=8))

#2. Feuchtigkeit
plot(density(df_weather$Luftfeuchtigkeit))
summary(df_weather$Luftfeuchtigkeit)

df_weather$Luftfeuchtigkeit_Groups <- as.numeric(cut2(df_weather$Luftfeuchtigkeit, g=8))

#3. Niederschlag
plot(density(df_weather$Niederschlag))
summary(df_weather$Niederschlag)

df_weather$Niederschlag_Groups <- as.numeric(cut2(df_weather$Niederschlag, g=4))

df_weather <- df_weather %>%
  select(Datum, Temperatur_Groups, MaxTemperatur_Groups, MinTemperatur_Groups, Luftfeuchtigkeit_Groups, Niederschlag_Groups)

df_weather$Datum <- as.Date(df_weather$Datum)

#An den Datensatz joinen
df <- left_join(df, df_weather, by ="Datum")

Machine Learning Features - Lag Variablen

filialen <- as.data.frame(unique(df$KundenNr))
df$Datum <- as.Date(df$Datum)

#Monat, Jahr, Tag, Quartal
df$day <- weekdays(as.Date(df$Datum))
df$monat <- month(as.Date(df$Datum))
df$Jahr <- year(as.Date(df$Datum))
df$quartal <- quarter(df$Datum, with_year = TRUE, fiscal_start = 1)


library(zoo)

test2 = df %>%
  group_by(ArtNr, KundenNr) %>%
  arrange(ArtNr, KundenNr , Datum) 


#Durchschnitt der letzten 2-7 Tage, 14 Tage und 30 Tage
test2 = test2 %>%
  mutate(VerkaufteMenge.lag1 = lag(VerkaufteMenge, n = 1)) %>%
  mutate(VerkaufteMenge.3.previous = rollapply(data = VerkaufteMenge.lag1, 
                                     width = 3, 
                                     FUN = mean, 
                                     align = "right", 
                                     fill = NA, 
                                     na.rm = T))
test2 = test2 %>%
  mutate(VerkaufteMenge.2.previous = rollapply(data = VerkaufteMenge.lag1, 
                                               width = 2, 
                                               FUN = mean, 
                                               align = "right", 
                                               fill = NA, 
                                               na.rm = T))
test2 = test2 %>%
  mutate(VerkaufteMenge.4.previous = rollapply(data = VerkaufteMenge.lag1, 
                                               width = 4, 
                                               FUN = mean, 
                                               align = "right", 
                                               fill = NA, 
                                               na.rm = T))

test2 = test2 %>%
  mutate(VerkaufteMenge.5.previous = rollapply(data = VerkaufteMenge.lag1, 
                                               width = 5, 
                                               FUN = mean, 
                                               align = "right", 
                                               fill = NA, 
                                               na.rm = T))

test2 = test2 %>%
  mutate(VerkaufteMenge.6.previous = rollapply(data = VerkaufteMenge.lag1, 
                                               width = 6, 
                                               FUN = mean, 
                                               align = "right", 
                                               fill = NA, 
                                               na.rm = T))
test2 = test2 %>%
  mutate(VerkaufteMenge.7.previous = rollapply(data = VerkaufteMenge.lag1, 
                                               width = 7, 
                                               FUN = mean, 
                                               align = "right", 
                                               fill = NA, 
                                               na.rm = T))

test2 = test2 %>%
  mutate(VerkaufteMenge.14.previous = rollapply(data = VerkaufteMenge.lag1, 
                                               width = 14, 
                                               FUN = mean, 
                                               align = "right", 
                                               fill = NA, 
                                               na.rm = T))

test2 = test2 %>%
  mutate(VerkaufteMenge.30.previous = rollapply(data = VerkaufteMenge.lag1, 
                                               width = 30, 
                                               FUN = mean, 
                                               align = "right", 
                                               fill = NA, 
                                               na.rm = T))




#Wochentage
test2$wochentag <- test2$day
test2$wochentag <- gsub("Montag",1, test2$wochentag)
test2$wochentag <- gsub("Dienstag",2, test2$wochentag)
test2$wochentag <- gsub("Mittwoch",3, test2$wochentag)
test2$wochentag <- gsub("Donnerstag",4, test2$wochentag)
test2$wochentag <- gsub("Freitag",5, test2$wochentag)
test2$wochentag <- gsub("Samstag",6, test2$wochentag)
test2$wochentag <- gsub("Sonntag",7, test2$wochentag)

weekend <- c(6,7)
test2$isWeekend <- if_else(test2$wochentag %in% weekend, 1,0)



#Vorwoche
test2 = test2 %>%
  group_by(ArtNr, KundenNr, wochentag) %>%
  arrange(wochentag, ArtNr, KundenNr , Datum) %>%
  mutate(Vorworche = lag(VerkaufteMenge, n = 1))
test2 = test2 %>%
  arrange(ArtNr, KundenNr , Datum)
#Verkaufszahlen der Vorwoche

#Vor 2 Wochen
test2 = test2 %>%
  group_by(ArtNr, KundenNr, wochentag) %>%
  arrange(wochentag, ArtNr, KundenNr , Datum) %>%
  mutate(VorVorworche = lag(VerkaufteMenge, n = 2))
test2 = test2 %>%
  arrange(ArtNr, KundenNr , Datum)

#Vor 3 Wochen
test2 = test2 %>%
  group_by(ArtNr, KundenNr, wochentag) %>%
  arrange(wochentag, ArtNr, KundenNr , Datum) %>%
  mutate(VorVorVorworche = lag(VerkaufteMenge, n = 3))
test2 = test2 %>%
  arrange(ArtNr, KundenNr , Datum)

#Vor 4 Wochen
test2 = test2 %>%
  group_by(ArtNr, KundenNr, wochentag) %>%
  arrange(wochentag, ArtNr, KundenNr , Datum) %>%
  mutate(VorVorVorVorworche = lag(VerkaufteMenge, n = 4))
test2 = test2 %>%
  arrange(ArtNr, KundenNr , Datum)

#Durchschnitt der Vorwochenwerte
test2$VorwochenAverage = rowMeans(test2[,c("Vorworche","VorVorworche", "VorVorVorworche", "VorVorVorVorworche")], na.rm=TRUE)
#Abweichung der Vorwoche zum Vorwochendurchschnitt
test2$VorwocheDifAverage = test2$Vorworche - test2$VorwochenAverage



test2 = test2 %>%
  group_by(ArtNr, KundenNr, isWeekend) %>%
  arrange(isWeekend, ArtNr, KundenNr , Datum) %>%
  mutate(WeekendAverage = (lag(Vorworche, n = 1) + lag(VorVorworche, n = 1) + lag(VorVorVorworche, n = 1) + lag(VorVorVorVorworche, n = 1) + lag(Vorworche, n = 2) + lag(VorVorworche, n = 2) + lag(VorVorVorworche, n = 2) + lag(VorVorVorVorworche, n = 2))/8)
test2 = test2 %>%
  arrange(ArtNr, KundenNr , Datum)
#dieses Feature nimmt immer den Durchschnitt der letzten 2 Einheiten. Also für Samstag = VorSonntag + VorSamstag, Für Sonntag = Gestern + VorSonntag
test2$WeekendAverage <- if_else(test2$isWeekend == 1, test2$WeekendAverage, 0)
#Dieses Feature steht nur für Wochenenden

test2 = test2 %>%
  group_by(ArtNr, KundenNr) %>%
  arrange(ArtNr, KundenNr , Datum) 
test2 = test2 %>%
  mutate(WeekendAverage2 = lag(WeekendAverage, n = 1))
test2$PreviousWeekendAverage <- if_else(test2$wochentag == 7, test2$WeekendAverage2, test2$WeekendAverage)
test2 <- subset(test2, select = - c(WeekendAverage2))



test2 = test2 %>%
  group_by(ArtNr, KundenNr, isWeekend) %>%
  arrange(isWeekend, ArtNr, KundenNr , Datum) %>%
  mutate(UnterwocheAverage = (lag(VorwochenAverage, n = 1) + lag(VorwochenAverage, n = 2) + lag(VorwochenAverage, n = 3) + lag(VorwochenAverage, n = 4) + lag(VorwochenAverage, n = 5))/5)
test2 = test2 %>%
  arrange(ArtNr, KundenNr , Datum)
#dieses Feature nimmt immer den Durchschnitt der letzten 5 Einheiten. Also für Freitag = Donnerstag + Mittwoch + Dienstag + Montag + Freitag
test2$UnterwocheAverage <- if_else(test2$isWeekend == 0, test2$UnterwocheAverage, 0)
#Dieses Feature steht nur für Unter der Woche

Feature Korrelation

library(corrplot)
library(caret)
df_cor <- as.data.frame(lapply(test2, as.numeric))
df_cor <- subset(df_cor, select = - c(Datum, day, ArtNr, KundenNr, Artikelgruppe))
df_cor <- na.omit(df_cor)
corMatrix <- cor(df_cor)
corrplot(corMatrix, type = "upper", order = "hclust", 
         tl.col = "black", tl.srt = 45)

One Hot Encoding

df <- test2

df$day <- as.factor(df$day)
df$ArtNr <- as.factor(df$ArtNr)
df$KundenNr <- as.factor(df$KundenNr)
df$Jahr <- as.factor(df$Jahr)
df$monat <- as.factor(df$monat)
df$quartal <- as.factor(df$quartal)

df$Temperatur_Groups <- as.factor(df$Temperatur_Groups)
df$MaxTemperatur_Groups <- as.factor(df$MaxTemperatur_Groups)
df$MinTemperatur_Groups <- as.factor(df$MinTemperatur_Groups)
df$Luftfeuchtigkeit_Groups <- as.factor(df$Luftfeuchtigkeit_Groups)
df$Niederschlag_Groups <- as.factor(df$Niederschlag_Groups)

df$holiday <- as.factor(df$holiday)


#Idee: Hier könnte man sowas wie A, B oder C Gruppe machen
df_tb <- setDT(df)
df_tb <- one_hot(df_tb, cols = "day")
df_tb <- one_hot(df_tb, cols = "ArtNr")
df_tb <- one_hot(df_tb, cols = "KundenNr")
df_tb <- one_hot(df_tb, cols = "Jahr")
df_tb <- one_hot(df_tb, cols = "monat")
df_tb <- one_hot(df_tb, cols = "quartal")

df_tb <- one_hot(df_tb, cols = "Temperatur_Groups")
df_tb <- one_hot(df_tb, cols = "MaxTemperatur_Groups")
df_tb <- one_hot(df_tb, cols = "MinTemperatur_Groups")
df_tb <- one_hot(df_tb, cols = "Luftfeuchtigkeit_Groups")
df_tb <- one_hot(df_tb, cols = "Niederschlag_Groups")

df_tb <- one_hot(df_tb, cols = "holiday")

df <- as.data.frame(df_tb)

Bereinigung

Bereinigung um NA’S

# Für die Lag variablen wird bei NA mit Median ersetzt

for(i in 87:ncol(df )){
  df [is.na(df [,i]), i] <- median(df [,i], na.rm = TRUE)
}

for(i in 9:16){
  df [is.na(df [,i]), i] <- 0
}

#NA Test
colnames(df)[colSums(is.na(df)) > 0]

for(i in 7:7){
  df [is.na(df [,i]), i] <- median(df [,i], na.rm = TRUE)
}

Modellierung der Nachfrage

Prophet Parameter Tuning

Für 4 Parameter werden hier die optimalen Parameter bestimmt. Idee: Mit jedem Parameter ein Modell fitten und anschließend den RMSE speichern. Nehme die Parameterkonstellation, welche den niedrigsten RMSE hat.

#einzelne Produkte
df_prophet <- df_final %>%
  filter(artnr == 1528) %>%
  select(datum, FinalerDemand)
colnames(df_prophet) <- c("ds","y")

#Logarithmieren
df_prophet$y <- log(df_prophet$y)
df_prophet <- df_prophet %>%
  filter(y != -Inf)

#Initial Trainingsset
l <- nrow(df_prophet) - 60
train_mc <- df_prophet[(1:l),]


#Hyper Parameter
cps <- 0.25
hpc <- 0.25
ys <- 10
ms <- 5

results <- data.frame(matrix(ncol = 5, nrow = 0))
colnames(results) <- c("RMSE","cps","hpc","ys","ms")


while (cps <= 1) {
  hpc <- 0.25
  while(hpc <= 1) {
    ys <- 10
    while (ys <= 25) {
      ms <- 5
      while (ms <= 10) {
     
    
    m = prophet(changepoint.prior.scale = cps, weekly.seasonality = TRUE, holidays = holidays,         holidays.prior.scale = hpc, yearly.seasonality = ys)
    m <- add_seasonality(m, name='monthly', period=30.5, fourier.order= ms)
    m <- fit.prophet(m, train_mc)
    
    future = as.data.frame(df_prophet$ds)
    colnames(future) <- c("ds")
    forecast = predict(m, future)
    
    val <- as.data.frame(forecast$yhat)
    val <- cbind(val, df_prophet$y)
    colnames(val) <- c("pred","actual")
    val$pred <- exp(val$pred)
    val$actual <- exp(val$actual)
    
    test <- val[-(1:l),]
    
    print(paste("RMSE Score: ",rmse(test$actual, test$pred), " ", cps, " ", hpc, sep = ""))
    
    result_schleife <- as.data.frame(rmse(test$actual, test$pred))
    colnames(result_schleife) <- c("RMSE")
    result_schleife$cps <- cps
    result_schleife$hpc <- hpc
    result_schleife$ys <- ys
    result_schleife$ms <- ms
    
    results <- rbind(results, result_schleife)
    
    
    ms <- ms + 2.5
  }
    
    ys <- ys + 5
  }
    
    
    hpc <- hpc + 0.25
  }
  
  cps <- cps + 0.25  
}

#Min Max rausbekommen
min <- min(results$RMSE)
max <- max(results$RMSE)

#Auf Minimum reduzieren
min_parameter <- results %>%
  filter(RMSE == min)
#Alternativer Ansatz: Sortieren und aus den besten den mit minimalen Paramtern wählen
newdata <- results[order(results$RMSE),]
newdata$rounded <- round(newdata$RMSE,1)

Cross-Validation

set.seed(305)
zeitpunkte_cv <- sample(seq(as.Date('2018-07-01'), as.Date('2019-07-01'), by="day"), 6)
zeitpunkte_cv <- sort(zeitpunkte_cv)
zeitpunkte_cv <- c(zeitpunkte_cv, as.Date("2019-07-14"))

library(Metrics)
library(MLmetrics)

Facebook Prophet Modellierung

for (value in zeitpunkte_cv) {

df_prophet <- df_final %>%
  filter(artnr == 1528) %>%
  select(datum, FinalerDemand)
colnames(df_prophet) <- c("ds","y")

df_prophet$y <- if_else(is.na(df_prophet$y) == TRUE, 0, df_prophet$y)
df_prophet <- df_prophet %>%
    filter(y > quantile(df_prophet$y, 0.01))


#Logarithmieren
df_prophet$y <- log(df_prophet$y)
df_prophet <- df_prophet %>%
  filter(y != -Inf)

cps <- 0.5
hpc <- 0.25
ys <- 10
ms <- 5
    
  
results <- data.frame(matrix(ncol = 3, nrow = 0))
colnames(results) <- c("Zeitpunkt","Actuals","Pred")

df_prophet <- df_prophet %>%
  filter(ds <= value)


i <- 0
while(i < 14) {
  l <- nrow(df_prophet) - 14 + i
  train_mc <- df_prophet[(1:l),]
  m = prophet(changepoint.prior.scale = cps, weekly.seasonality = TRUE, holidays = holidays,         holidays.prior.scale = hpc, yearly.seasonality = ys)
    m <- add_seasonality(m, name='monthly', period=30.5, fourier.order= ms)
    m <- fit.prophet(m, train_mc)
    future = as.data.frame(df_prophet$ds)
    future = as.data.frame(future[(l:l+1),])
    colnames(future) <- c("ds")
    forecast = predict(m, future)
    
    val <- as.data.frame(forecast$yhat)
    actual = as.data.frame(df_prophet[(l:l+1),])
    val$actual <- actual$y
    colnames(val) <- c("pred","actual")
    val$pred <- exp(val$pred)
    val$actual <- exp(val$actual)

    result_schleife <- as.data.frame(actual$ds)
    colnames(result_schleife) <- c("Zeitpunkt")
    result_schleife$Actuals <- val$actual
    result_schleife$Pred <- val$pred

    results <- rbind(results, result_schleife)
    
    
    
    
  
  i <- i +1
}

print(rmse(results$Actuals, results$Pred))
print(MAPE(results$Actuals, results$Pred))
print(MAE(results$Actuals, results$Pred))

}    

XGBoost Modellierung

df <- df[order(as.Date(df$Datum)),]


library(data.table)
library(mlr)
library(xgboost)

results <- data.frame(matrix(ncol = 4, nrow = 0))
colnames(results) <- c("Zeitpunkt","RMSE","MAE","MAPE")



for (value in zeitpunkte_cv) {
  
  df_xgb <- df %>%
    filter(Datum <= as.Date(value)) %>%
    filter(ArtNr_1525 == 1)
  
  df_xgb <- df_xgb %>%
    filter(VerkaufteMenge > quantile(df_xgb$VerkaufteMenge, 0.01))
  #Die ganz krassen Ausreißer sollen nicht enthalten sein
  
  l <- nrow(df_xgb) - 14
  #train_set <- df_xgb[(1:l),]
  train_set <- df_xgb[(1:(l-30)),]
  test_set <- df_xgb[((l-29):l),]
  val_set <- df_xgb[((l+1):nrow(df_xgb)),]
  
  train_x <- setDT(train_set)
  test_x <- setDT(test_set)
  val_x <- setDT(val_set)
  
  labels_train <- train_x$VerkaufteMenge
  labels_test <- test_x$VerkaufteMenge
  labels_val <- val_x$VerkaufteMenge

  train_x <- subset(train_x, select = - c(VerkaufteMenge, GelieferteMenge, Datum, Artikelgruppe))
  test_x <- subset(test_x, select = - c(VerkaufteMenge, GelieferteMenge, Datum, Artikelgruppe))
  val_x <- subset(val_x, select = - c(VerkaufteMenge, GelieferteMenge, Datum, Artikelgruppe))
  
  train_x <- data.matrix(train_x, rownames.force = NA)
  test_x <- data.matrix(test_x, rownames.force = NA)
  val_x <- data.matrix(val_x, rownames.force = NA)
  
  dtrain <- xgb.DMatrix(data = train_x,label = labels_train) 
  dtest <- xgb.DMatrix(data = test_x,label=labels_test)
  dval <- xgb.DMatrix(data = val_x,label=labels_val)
  
  watchlist <- list(train=dtrain, test=dtest)

  set.seed(1)
  st <- xgb.train(data=dtrain, booster = "gbtree", max.depth=7, lambda = 2, alpha = 2, colsample_bytree = 0.5, min_child_weight = 0.5, early_stopping_rounds = 200, nthread = 2, gamma = 0, eta= 0.01, nrounds=10000, watchlist=watchlist, eval.metric = "rmse", objective = "reg:linear")
  
  xgbpred <- predict (st,dval)
  
  result_schleife <- as.data.frame(as.Date(value))
  colnames(result_schleife) <- c("Zeitpunkt")
  result_schleife$RMSE <- rmse(labels_val, round(xgbpred,0))
  result_schleife$MAE <- MAE(round(xgbpred,0), labels_val)
  result_schleife$MAPE <- MAPE(round(xgbpred,0), labels_val)

    
  results <- rbind(results, result_schleife)
  

}


mean(results$RMSE)
mean(results$MAE)
mean(results$MAPE)

Neural Net Modellierung

library(keras)

results <- data.frame(matrix(ncol = 4, nrow = 0))
colnames(results) <- c("Zeitpunkt","RMSE","MAE","MAPE")


for (value in zeitpunkte_cv) {
  
  df_nn <- df %>%
    filter(Datum <= as.Date(value)) %>%
    filter(ArtNr_1526 == 1)
  
   df_nn <- df_nn %>%
    filter(VerkaufteMenge > quantile(df_nn$VerkaufteMenge, 0.01))
  #Die ganz krassen Ausreißer sollen nicht enthalten sein

  l <- nrow(df_nn) - 14
  
  labels <- df_nn$VerkaufteMenge
  
  df_nn <- subset(df_nn, select = - c(VerkaufteMenge, GelieferteMenge, Artikelgruppe, Datum))
  
  df_nn[,(ncol(df_nn) - 8):(ncol(df_nn))] <- scale(df_nn[,(ncol(df_nn) - 8):(ncol(df_nn))])
  df_nn[,(ncol(df_nn) - 19):(ncol(df_nn) - 11)] <- scale(df_nn[,(ncol(df_nn) - 19):(ncol(df_nn) - 11)])
  
  
  train_data <- df_nn[(1:(l-30)),]
  test_data <- df_nn[((l-29):l),]
  val_data <- df_nn[((l+1):nrow(df_nn)),]

  train_labels <- labels[(1:(l-30))]
  test_labels <- labels[((l-29):l)]
  val_labels <- labels[((l+1):length(labels))]
  
  train_data  <- as.data.frame(lapply(train_data, as.numeric))
  test_data  <- as.data.frame(lapply(test_data, as.numeric))
  val_data  <- as.data.frame(lapply(val_data, as.numeric))

  train_data <- as.matrix(train_data)
  test_data <- as.matrix(test_data)
  val_data <- as.matrix(val_data)
  
  
  build_model <- function() {
  
  model <- keras_model_sequential() %>%
    layer_dense(units = 64, activation = "relu", input_shape = dim(train_data)[2]) %>%
    layer_dropout(rate=0.4) %>%
    layer_dense(units = 32, activation = "relu") %>%
    layer_dropout(rate=0.2) %>%
    layer_dense(units = 16, activation = "relu") %>%
    layer_dense(units = 1)
  
  model %>% compile(
    loss = 'mean_squared_error',
    optimizer = optimizer_rmsprop(),
    metrics = list("mean_absolute_error")
  )
  
  model
  }

  model <- build_model()
  model %>% summary()
  
  
  print_dot_callback <- callback_lambda(
    on_epoch_end = function(epoch, logs) {
      if (epoch %% 80 == 0) cat("\n")
      cat(".")
    }
  )    
  
  epochs <- 500
 
  early_stop <- callback_early_stopping(monitor = "val_loss", patience = 20)
  
  set.seed(41)
  
  model <- build_model()
  history <- model %>% fit(
    train_data,
    train_labels,
    epochs = epochs,
    validation_split = 0.2,
    verbose = 0,
    batch_size = 32,
    callbacks = list(early_stop, print_dot_callback)
  )
  
  test_predictions <- model %>% predict(test_data)

  result_schleife <- as.data.frame(as.Date(value))
  colnames(result_schleife) <- c("Zeitpunkt")
  result_schleife$RMSE <- rmse(test_labels, round(test_predictions[ , 1],0))
  result_schleife$MAE <- MAE(round(test_predictions[ , 1],0), test_labels)
  result_schleife$MAPE <- MAPE(round(test_predictions[ , 1],0), test_labels)

    
  results <- rbind(results, result_schleife)
    
  
  
}  

Residual Analyse - Prophet

df_prophet <- df_final %>%
  filter(artnr == 1526) %>%
  select(datum, FinalerDemand)
colnames(df_prophet) <- c("ds","y")

df_prophet$y <- if_else(is.na(df_prophet$y) == TRUE, 0, df_prophet$y)
df_prophet <- df_prophet %>%
    filter(y > quantile(df_prophet$y, 0.01))


#Logarithmieren
df_prophet$y <- log(df_prophet$y)
df_prophet <- df_prophet %>%
  filter(y != -Inf)

cps <- 0.5
hpc <- 0.25
ys <- 10
ms <- 5
    
  
results <- data.frame(matrix(ncol = 3, nrow = 0))
colnames(results) <- c("Zeitpunkt","Actuals","Pred")

df_prophet <- df_prophet %>%
  filter(ds <= as.Date("2019-06-30"))


i <- 0
while(i < 60) {
  l <- nrow(df_prophet) - 60 + i
  train_mc <- df_prophet[(1:l),]
  m = prophet(changepoint.prior.scale = cps, weekly.seasonality = TRUE, holidays = holidays,         holidays.prior.scale = hpc, yearly.seasonality = ys)
    m <- add_seasonality(m, name='monthly', period=30.5, fourier.order= ms)
    m <- fit.prophet(m, train_mc)
    future = as.data.frame(df_prophet$ds)
    future = as.data.frame(future[(l:l+1),])
    colnames(future) <- c("ds")
    forecast = predict(m, future)
    
    val <- as.data.frame(forecast$yhat)
    actual = as.data.frame(df_prophet[(l:l+1),])
    val$actual <- actual$y
    colnames(val) <- c("pred","actual")
    val$pred <- exp(val$pred)
    val$actual <- exp(val$actual)

    result_schleife <- as.data.frame(actual$ds)
    colnames(result_schleife) <- c("Zeitpunkt")
    result_schleife$Actuals <- val$actual
    result_schleife$Pred <- val$pred

    results <- rbind(results, result_schleife)
    
    
    
    
  
  i <- i +1
}

print(rmse(results$Actuals, results$Pred))
print(MAPE(results$Actuals, results$Pred))
print(MAE(results$Actuals, results$Pred))


#
#Residuals
results$residuals <- results$Pred - results$Actuals
#Mean von von den Fehlern
mean(results$residuals)
#Prediction um mean anpassen
results$predAdjusted <- results$Pred - mean(results$residuals)
#Jetzt liegt der mean bei 0
round(mean(results$predAdjusted - results$Actuals),1)
#Finale Fehler
results$residualsAdjusted <- results$predAdjusted - results$Actuals

#Standardabweichung der Fehler
sd_prophet <- sd(results$residualsAdjusted)

#Fehler Normalverteilung
set.seed(601)
y <- rnorm(100000, mean = 0, sd = sd_prophet)
hist(y, breaks = 100)

y <- round(y,0)

y <- as.data.frame(y)
colnames(y) <- c("y")
residaul_prophet <- count(y,y)
residaul_prophet$prob <- residaul_prophet$n / sum(residaul_prophet$n)
sum(residaul_prophet$prob)

Prophet Scenario Generation

zeitpunkte_cv <- c(as.Date("2019-07-14"))

for (value in zeitpunkte_cv) {

df_prophet <- df_final %>%
  filter(artnr == 1526) %>%
  select(datum, FinalerDemand)
colnames(df_prophet) <- c("ds","y")

df_prophet$y <- if_else(is.na(df_prophet$y) == TRUE, 0, df_prophet$y)
df_prophet <- df_prophet %>%
    filter(y > quantile(df_prophet$y, 0.01))


#Logarithmieren
df_prophet$y <- log(df_prophet$y)
df_prophet <- df_prophet %>%
  filter(y != -Inf)

cps <- 0.25
hpc <- 0.25
ys <- 10
ms <- 5
    
  
results <- data.frame(matrix(ncol = 3, nrow = 0))
colnames(results) <- c("Zeitpunkt","Actuals","Pred")

df_prophet <- df_prophet %>%
  filter(ds <= value)


i <- 0
while(i < 14) {
  l <- nrow(df_prophet) - 14 + i
  train_mc <- df_prophet[(1:l),]
  m = prophet(changepoint.prior.scale = cps, weekly.seasonality = TRUE, holidays = holidays,         holidays.prior.scale = hpc, yearly.seasonality = ys)
    m <- add_seasonality(m, name='monthly', period=30.5, fourier.order= ms)
    m <- fit.prophet(m, train_mc)
    future = as.data.frame(df_prophet$ds)
    future = as.data.frame(future[(l:l+1),])
    colnames(future) <- c("ds")
    forecast = predict(m, future)
    
    val <- as.data.frame(forecast$yhat)
    actual = as.data.frame(df_prophet[(l:l+1),])
    val$actual <- actual$y
    colnames(val) <- c("pred","actual")
    val$pred <- exp(val$pred)
    val$actual <- exp(val$actual)

    result_schleife <- as.data.frame(actual$ds)
    colnames(result_schleife) <- c("Zeitpunkt")
    result_schleife$Actuals <- val$actual
    result_schleife$Pred <- val$pred

    results <- rbind(results, result_schleife)
    
    
    
    
  
  i <- i +1
}

print(rmse(results$Actuals, results$Pred))
print(MAPE(results$Actuals, results$Pred))
print(MAE(results$Actuals, results$Pred))



}

predictions_prophet <- round(results$Pred,0)
predictions_colnames <- results$Zeitpunkt

for (value in predictions_prophet) {

residaul_prophet <- cbind(residaul_prophet, value + residaul_prophet$y)  
   
  
}  

colnames(residaul_prophet) <- c("y", "n", "prob", c(predictions_colnames))


residaul_prophet$index <- seq(1, nrow(residaul_prophet), by = 1)
residaul_prophet$Sc <- paste("S", residaul_prophet$index, sep = "")
write.csv(residaul_prophet, "C:/Users/Denni/Documents/R 3.5 Files/Masterarbeit/Untitled Folder/Prophet_Laugenbroetchen.csv")

Residual Analyse - XGBoost

df <- df[order(as.Date(df$Datum)),]

library(data.table)
library(mlr)
library(xgboost)

results <- data.frame(matrix(ncol = 3, nrow = 0))
colnames(results) <- c("Zeitpunkt","Actuals","Pred")


for (value in zeitpunkte_cv) {
  
  df_xgb <- df %>%
    filter(Datum <= as.Date(value)) %>%
    filter(ArtNr_1526 == 1)
  
  df_xgb <- df_xgb %>%
    filter(VerkaufteMenge > quantile(df_xgb$VerkaufteMenge, 0.01))
  #Die ganz krassen Ausreißer sollen nicht enthalten sein
  
  l <- nrow(df_xgb) - 14
  #train_set <- df_xgb[(1:l),]
  train_set <- df_xgb[(1:(l-30)),]
  test_set <- df_xgb[((l-29):l),]
  val_set <- df_xgb[((l+1):nrow(df_xgb)),]
  
  train_x <- setDT(train_set)
  test_x <- setDT(test_set)
  val_x <- setDT(val_set)
  
  labels_train <- train_x$VerkaufteMenge
  labels_test <- test_x$VerkaufteMenge
  labels_val <- val_x$VerkaufteMenge

  train_x <- subset(train_x, select = - c(VerkaufteMenge, GelieferteMenge, Datum, Artikelgruppe))
  test_x <- subset(test_x, select = - c(VerkaufteMenge, GelieferteMenge, Datum, Artikelgruppe))
  val_x <- subset(val_x, select = - c(VerkaufteMenge, GelieferteMenge, Datum, Artikelgruppe))
  
  train_x <- data.matrix(train_x, rownames.force = NA)
  test_x <- data.matrix(test_x, rownames.force = NA)
  val_x <- data.matrix(val_x, rownames.force = NA)
  
  dtrain <- xgb.DMatrix(data = train_x,label = labels_train) 
  dtest <- xgb.DMatrix(data = test_x,label=labels_test)
  dval <- xgb.DMatrix(data = val_x,label=labels_val)
  
  watchlist <- list(train=dtrain, test=dtest)

  set.seed(1)
  st <- xgb.train(data=dtrain, booster = "gbtree", max.depth=7, lambda = 2, alpha = 2, colsample_bytree = 0.5, min_child_weight = 0.5, early_stopping_rounds = 200, nthread = 2, gamma = 0, eta= 0.01, nrounds=10000, watchlist=watchlist, eval.metric = "rmse", objective = "reg:linear")
  
  xgbpred <- predict(st,dtest)
  
  
  

}
results <- as.data.frame(test_set$Datum)
results$Actuals <- labels_test
results$Pred <- xgbpred
colnames(results) <- c("Zeitpunkte", "Actuals", "Pred")
rmse(labels_test,xgbpred)


#Residuals
results$residuals <- results$Pred - results$Actuals
#Mean von von den Fehlern
mean(results$residuals)
#Prediction um mean anpassen
results$predAdjusted <- results$Pred - mean(results$residuals)
#Jetzt liegt der mean bei 0
round(mean(results$predAdjusted - results$Actuals),1)
#Finale Fehler
results$residualsAdjusted <- results$predAdjusted - results$Actuals

#Standardabweichung der Fehler
sd_xgb <- sd(results$residualsAdjusted)


#Fehler Normalverteilung
set.seed(602)
y <- rnorm(100000, mean = 0, sd = sd_xgb)
hist(y, breaks = 100)

y <- round(y,0)

y <- as.data.frame(y)
colnames(y) <- c("y")
residaul_xgb <- count(y,y)
residaul_xgb$prob <- residaul_xgb$n / sum(residaul_xgb$n)
sum(residaul_xgb$prob)

Scenario Generation - XGB

xgbpred <- predict(st,dval)

results <- as.data.frame(val_set$Datum)
results$Actuals <- labels_val
results$Pred <- xgbpred
colnames(results) <- c("Zeitpunkt", "Actuals", "Pred")
rmse(labels_val,xgbpred)




predictions_xgb <- round(results$Pred,0)
predictions_colnames <- results$Zeitpunkt

for (value in predictions_xgb) {

residaul_xgb <- cbind(residaul_xgb, value + residaul_xgb$y)  
   
  
}  

colnames(residaul_xgb) <- c("y", "n", "prob", c(predictions_colnames))


residaul_xgb$index <- seq(1, nrow(residaul_xgb), by = 1)
residaul_xgb$Sc <- paste("S", residaul_xgb$index, sep = "")
write.csv(residaul_xgb, "C:/Users/Denni/Documents/R 3.5 Files/Masterarbeit/Untitled Folder/xgb_Laugenbroetchen.csv")

Residual Analyse - Neural Net

library(keras)

results <- data.frame(matrix(ncol = 3, nrow = 0))
colnames(results) <- c("Zeitpunkt","Actuals","Pred")


for (value in zeitpunkte_cv) {
  
  df_nn <- df %>%
    filter(Datum <= as.Date(value)) %>%
    filter(ArtNr_1526 == 1)
  
   df_nn <- df_nn %>%
    filter(VerkaufteMenge > quantile(df_nn$VerkaufteMenge, 0.01))
  #Die ganz krassen Ausreißer sollen nicht enthalten sein

  l <- nrow(df_nn) - 14
  
  labels <- df_nn$VerkaufteMenge
  
  df_nn <- subset(df_nn, select = - c(VerkaufteMenge, GelieferteMenge, Artikelgruppe, Datum))
  
  df_nn[,(ncol(df_nn) - 8):(ncol(df_nn))] <- scale(df_nn[,(ncol(df_nn) - 8):(ncol(df_nn))])
  df_nn[,(ncol(df_nn) - 19):(ncol(df_nn) - 11)] <- scale(df_nn[,(ncol(df_nn) - 19):(ncol(df_nn) - 11)])
  
  
  train_data <- df_nn[(1:(l-30)),]
  test_data <- df_nn[((l-29):l),]
  val_data <- df_nn[((l+1):nrow(df_nn)),]

  train_labels <- labels[(1:(l-30))]
  test_labels <- labels[((l-29):l)]
  val_labels <- labels[((l+1):length(labels))]
  
  train_data  <- as.data.frame(lapply(train_data, as.numeric))
  test_data  <- as.data.frame(lapply(test_data, as.numeric))
  val_data  <- as.data.frame(lapply(val_data, as.numeric))

  train_data <- as.matrix(train_data)
  test_data <- as.matrix(test_data)
  val_data <- as.matrix(val_data)
  
  
  build_model <- function() {
  
  model <- keras_model_sequential() %>%
    layer_dense(units = 64, activation = "relu", input_shape = dim(train_data)[2]) %>%
    layer_dropout(rate=0.4) %>%
    layer_dense(units = 32, activation = "relu") %>%
    layer_dropout(rate=0.2) %>%
    layer_dense(units = 16, activation = "relu") %>%
    layer_dense(units = 1)
  
  model %>% compile(
    loss = 'mean_squared_error',
    optimizer = optimizer_rmsprop(),
    metrics = list("mean_absolute_error")
  )
  
  model
  }

  model <- build_model()
  model %>% summary()
  
  
  # Display training progress by printing a single dot for each completed epoch.
  print_dot_callback <- callback_lambda(
    on_epoch_end = function(epoch, logs) {
      if (epoch %% 80 == 0) cat("\n")
      cat(".")
    }
  )    
  
  epochs <- 500
 
  early_stop <- callback_early_stopping(monitor = "val_loss", patience = 20)
  
  set.seed(41)
  
  model <- build_model()
  history <- model %>% fit(
    train_data,
    train_labels,
    epochs = epochs,
    validation_split = 0.2,
    verbose = 0,
    batch_size = 32,
    callbacks = list(early_stop, print_dot_callback)
  )
  
  test_predictions <- model %>% predict(test_data)

  
}

results <- as.data.frame(seq(1, 30, by = 1))
results$Actuals <- test_labels
results$Pred <- test_predictions[ , 1]
colnames(results) <- c("Zeitpunkte", "Actuals", "Pred")
rmse(test_labels,test_predictions[ , 1])


#Residuals
results$residuals <- results$Pred - results$Actuals
#Mean von von den Fehlern
mean(results$residuals)
#Prediction um mean anpassen
results$predAdjusted <- results$Pred - mean(results$residuals)
#Jetzt liegt der mean bei 0
round(mean(results$predAdjusted - results$Actuals),1)
#Finale Fehler
results$residualsAdjusted <- results$predAdjusted - results$Actuals

#Standardabweichung der Fehler
sd_nn <- sd(results$residualsAdjusted)

#Normalverteilung der Fehler
set.seed(604)
y <- rnorm(100000, mean = 0, sd = sd_nn)
hist(y, breaks = 100)

y <- round(y,0)

y <- as.data.frame(y)
colnames(y) <- c("y")
residaul_nn <- count(y,y)
residaul_nn$prob <- residaul_nn$n / sum(residaul_nn$n)
sum(residaul_nn$prob)

Scenario Generation - Neural Net

test_predictions <- model %>% predict(val_data)

results <- as.data.frame(seq(1, 14, by = 1))
results$Actuals <- val_labels
results$Pred <- test_predictions[ , 1]
colnames(results) <- c("Zeitpunkt", "Actuals", "Pred")
rmse(results$Actuals,results$Pred)




predictions_nn <- round(results$Pred,0)
predictions_colnames <- results$Zeitpunkt

for (value in predictions_nn) {

residaul_nn <- cbind(residaul_nn, value + residaul_nn$y)  
   
  
}  

colnames(residaul_nn) <- c("y", "n", "prob", c(predictions_colnames))


residaul_nn$index <- seq(1, nrow(residaul_nn), by = 1)
residaul_nn$Sc <- paste("S", residaul_nn$index, sep = "")
write.csv(residaul_nn, "C:/Users/Denni/Documents/R 3.5 Files/Masterarbeit/Untitled Folder/nn_Laugenbroetchen.csv")